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INTRODUCTION 


The purpose of this thesis is to compare the fourth-order Runge-Kutta 
numerical integration technique with two of Dr. E. R. Shanks 1 numerical 
integration formulas in a mathematical model used for computing rigid body, 
first-stage trajectories of a Saturn space vehicle. The scope of the problem is 
to describe the mathematical model and to present a derivation of the differen- 
tial equations of motion which it comprises; to establish a basis for comparing 
the integration techniques; and to generate sufficient comparison data to 
establish which integration technique is the most practical to use in this model 
from a standpoint of computer run time and accuracy. 

As an employee of the National Aeronautics and Space Administration, 
the author was required to develop a computer program for simulating the 
first-stage flight of a Saturn space vehicle. The resulting program requires 
that six first-order and three second-order differential equations be numeri- 
cally solved to predict the translational and rotational motion of a space 
vehicle. An efficient integration technique had to be chosen that would provide 
a solution within the accuracy of the data which describe the vehicle dynamic 
response characteristics . In a search for the most acceptable integration 
technique, the fourth -order Runge-Kutta formula and two of Dr. E. B. Shanks f 



integration formulas were each used in the program to compute the first-stage 
trajectory of a typical Saturn vehicle. Comparative data were generated for 
each formula over a range of integration step sizes. These data are analyzed 
in Chapter V. 

The first four chapters define the mathematical model in which the 
integration formulas are compared. This model consists of several smaller 
models that simulate the vehicle subsystems and the physical properties of the 
earth. These subsystems and the earth’s physical characteristics are 
described in Chapter I. The different coordinate systems to which the model 
refers motion and the related transformations are defined in Chapter II, and 
the differential equations of motion are derived in Chapter ILL In Chapter IV, 
a definition of the mathematical model is completed with a description of the 
forces and moments that act upon a space vehicle during flight. 
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CHAPTER I 


SATURN VEHICLE SUBSYSTEMS AND PHYSICAL PROPERTIES 
OF THE EARTH TO BE SIMULATED 

In order to design and build an extremely complex and expensive vehicle 
such as the Saturn V, a means of simulating or predicting how a space vehicle 
with specific subsystems and dynamic characteristics will react during actual 
flight is mandatory. To this end, mathematical models are developed for each 
subsystem and then combined into a total model of the vehicle. Models of the 
earth’s shape, gravity field, and atmosphere are merged with the vehicle 
model to establish a two-body system. The total model is then programmed on 
an electronic computer and used for predicting the response of a vehicle to 
various atmospheric conditions and vehicle tolerances. 

The Saturn V subsystems that guide and control the vehicle will be 
described first. A knowledge of these subsystems and their functions is 
necessary before a model can be developed to simulate them. 

In this paper the terms navigation, guidance, and control are defined 
as follows: 

Navigation is the determination of the vehicle’s present position and 

velocity from measurements made onboard the vehicle. 
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Guidance is the computations of maneuvers necessary to achieve the 
desired flight path. 

Control is the execution of the guidance maneuver by controlling the 
direction of part of the vehicles thrust. 

The Saturn V launch vehicle is guided during powered flight by naviga- 
tion, guidance, and control equipment located in the Instrument Unit (IU) . An 
inertial platform is used to provide three orthogonal space -stabilized directions 
for acceleration and attitude measurements. A launch vehicle digital computer 
(LVDC) solves guidance equations, and an analog flight control computer 
executes flight control functions. Analog signals are converted to digital 
numbers and vice versa by means of a launch vehicle data adapter (LVDA) . 

The vehicle -fixed coordinate system and the Saturn sign convention for 
control variables are presented in Figure 1. Rotational motion about the 

axis is called roll motion, rotational motion about the Y axis is called 

m 

pitch motion* and rotational motion about the axis is referred to as yaw 

motion. 

Figure 2 presents the Saturn V inertial platform configuration. The 
gimbal system allows the vehicle to rotate freely without disturbing the gyro- 
stabilized inertial gimbal. An orthogonal, right-handed, space -fixed coordi- 
nate system (X, Y, Z) is established by the input axes of three single -degree - 
of -freedom gyroscopes. Acceleration and attitude measurements are taken 
with respect to this coordinate system. Three integrating accelerometers, 
mounted on the platform’s inertial gimbal, measure the three components of 
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velocity. Vehicle attitude is measured by resolvers located at the gimbal 
pivot points. 

The three angles measured by the inertial platform resolvers are 
classically known as Eulerian angles, a discussion of which should aid in 
understanding the function of the platform. 

Eulerian angles evolved from the problem of describing the orientation 
of a rigid body that is free to turn about a point 0. One solution would be to 
assign the coordinates of two particles in the body, not on the same line 
through 0. Because this method involves the use of six parameters that cannot 
be varied independently, it is rejected. A much simpler method is to define 
the orientation of the body with respect to a reference frame by the use of three 
Eulerian angles. 

Figure 3 shows two orthogonal right-handed triads (X, Y, Z) and 

(X , Y , Z ) with their origins at point 0. The triad (X , Y , Z ) is 
m m m' & m m m 

fixed in a rigid body that is free to turn about 0, and the triad (X, Y, Z) 

serves as a frame of reference. Assuming that the (X , Y , Z ) triad 

m m m 

initially coincides with the reference triad, it can be revolved to the orienta- 
tion shows in Figure 3 by the three following right-hand rotations. 

1. First rotate about Y through <p . This brings the movable triad 

(X , Y , Z ) into coincidence with (X ! , Y, Z T ) . 
m m m 

2. Next, rotate about Z ! through w . This brings (X , Y , Z ) 

y & m m m 

into coincidence with (X , Y\ Z ! ) . 

m 

3. Finally, rotate about X^ through <p . This brings 

(X , Y , Z ) into the final position, 
m m m 


7 




8 



I 

Note that all possible positions of the body can be obtained by assigning 
values to <p ^ 9 cp^, and cp^ in the ranges 0 < 0^ < 27r , 0 < c p < 7r , and 

0 < cp < 2tt . 

The symbols used in the discussion of Euler ian angles were chosen to 

coincide with those used in Figures i and 2. That is, the (X, Y, Z) triad is 

the one established by the vehicle inertial platform and the (X , Y , Z ) 

m m m 

triad corresponds to the vehicle-fixed triad in Figure i. An inspection of 

Figure 2 reveals that the Eulerian angles <p^ f cp^, cp ^ are precisely the 

angles measured by the inertial platform resolvers. 

Attitude control is achieved during the S-IC stage of powered flight by 

gimbaling or deflecting the stage’s four outboard engines (engines i through 4 

in Figure 1) . Each of these four engines has two hydraulic actuators for 

gimbaling the engine about its stage attach point. To determine the amount 

these engines are to be gimbaled, the LVDC compares the instantaneous 

vehicle attitude with the attitude computed by the guidance scheme. Attitude 

error signals (0 , 0 , ip^) are derived from the difference between the 

platform gimbal angles <p , cp and the desired attitude angles 

( X > X , X ) • The vehicle angular velocities about the X , Y , and Z 
VA p’ y’ 'Y & m m m 

axes are measured by three rate gyroscopes located in the IU. In the flight 

control computer, the signals from the rate gyroscopes ( $ ,0 , 0 ) and 

p y r 

the attitude error signals are filtered, multiplied by a gain factor, and then 
combined to generate the control commands for the engine actuators. The 
attitude error equations and the control laws for combining the filtered attitude 
errors and vehicle angular rates are presented in Chapter IV. 
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During the S-IC stage of flight, the vehicle must execute pitch, yaw, 
and roll maneuvers. A yaw maneuver during the first critical seconds of 
flight guides the vehicle away from the launch umbilical tower. The axis 

is aligned along the flight azimuth by the roll maneuver. The pitch maneuver 
is designed to guide the vehicle along a flight path that results in minimum 
aerodynamic structural loading. This maneuver allows the vehicle to build a 
component of velocity tangent to the earth's surface as it gains altitude. All 
three maneuvers are preprogrammed in the LVDC as time -dependent poly- 
nomials. 

The LVDC contains a closed-loop guidance scheme for providing 
guidance commands during the S-II and S-IVB stages of powered flight. This 
scheme determines attitude commands from the components of position and 
velocity of the vehicle center of gravity (CG) relative to a space-fixed coordi- 
nate system orientated as the (X, Y, Z) platform coordinate system. Initial 
velocity of the vehicle because of the earth's rotation and the platform velocity 
from the integrating accelerometers are algebraically summed by the LVDC 
to obtain space-fixed velocity of the vehicle. Position is determined by 
integrating the space -fixed vehicle velocity. 

With the description of the major vehicle subsystems complete, a 
presentation of the models of the earth's shape, motion, gravity field, and 
atmosphere can now be given. The standard model for the earth's shape is 
known as the Fischer Ellipsoid, a model which considers the earth to be an 
elliptical spheroid ( Figure 4) . The intersection of the earth and a plane that 
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THE FISCHER ELLIPSOID 


contains the earth T s polar axis is an ellipse; the intersection of the earth and 
a plane that is perpendicular to the polar axis is a circle. 

The only motion of the earth that will be considered is rotational; that 
is. the center of the earth will be assumed fixed in space with the earth 
rotating about its polar axis. Translation of the earth about the sun and 
precession of the earth’s polar axis during the year have no significant effect 
on the motion of a space vehicle with respect to the earth. 

A Newtonian potential function is used to predict the magnitude and 
direction of the earth’s gravity at any point above the earth’s surface. The 
oblate potential of the earth is assumed to contain the second, third, and 
fourth spherical harmonics. A discussion of this potential is lengthy and will 
not be presented here. A complete treatment of the subject is given in Chapters 
I and II of Reference 1. 

A model reference atmosphere for Cape Kennedy, Florida (Reference 
2) , which provides atmospheric information, is based on annual median values 
of measured atmospheric parameters. Multisegment polynomials define each 
atmospheric parameter as a function of altitude. 
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CHAPTER II 


REFERENCE COORDINATE SYSTEMS AND TRANSFORMATIONS 

In order to determine the motion of a vehicle with respect to the 
rotating earth, six reference three-dimensional coordinate systems are used. 
Each system is defined to be right-handed and orthogonal to simplify computa- 
tion. Rotational motion of a vehicle is computed in the vehicle coordinate 
system (X , Y , Z^) defined in Figure 2. The other five systems consist 
of two earth -fixed and three space -fixed systems. Here, the terms "space - 
fixed" and "earth-fixed" refer to coordinate systems that have constant orien- 
tation and position with respect to the solar system and the earth, respectively 
The three space-fixed systems (Figure 5) are referred to as the surface triad 

(X , Y , Z ) , the earth-centered triad (X , Y , Z ) , and the equatorial 
s s s c c c ^ 

triad (U, V, W) . The two earth-fixed systems (Figure 6) are referred to 

as the launch-point triad (X , Y , Z ) and the earth -equatorial triad 

e e e 

(U , V , W ) . 
x e’ e e 

At the instant the vehicle inertial platform is released, the origin of 

the (X g , Y g , Z ) triad is located at the launch point. As time passes, the 

surface of the earth moves beneath the (X , Y , Z ) origin because of the 

s s s 
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SPACE- FIXED COORDINATE SYSTEMS 




earth ! s rotation. This system is defined to have the same orientation as the 


(X, Y, Z) space-fixed triad established by the inertial platform (Figure 2) . 

At release, the inertial platform X axis is aligned to the local vertical; the 

Z axis is tangent to the earth ! s surface at the launch site and lies in the flight 

plane; and the Y axis forms a right-hand system. The (X , Y , Z ) triad 

s s s 

is defined in this manner to simulate the function of the vehicle inertial 
platform . 

The earth-centered triad (X , Y , Z ) has the same orientation as 
the surface triad; its origin is located at the center of the earth as the name 
implies . Some of the initial conditions and certain trajectory parameters are 
computed with respect to this system. 

A second earth -centered triad (U, V, W) , the equatorial system, is 
defined to simplify computations that involve the earth T s rotational velocity. 
The U axis, which is directed along the earth T s rotational velocity vector, 
passes through the North Pole and is perpendicular to the equatorial plane. 

At the instant the vehicle inertial platform is released, the W axis lies in the 
plane of the launch meridian. 

One earth -fixed triad (X , Y , Z ) is used by the model for comput- 

e e e 

ing motion of the vehicle with respect to the launch point ( Figure 6) . This 

triad has its origin at the launch site and is coincident with the (X , Y , Z ) 

s s s 

triad at the instant the vehicle platform is released. 

A second earth -fixed triad (U , V , W ) has its origin at the center 

e e © 

of the earth and is coincident with the (U, V, W) triad at platform release. 
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This system is used in deriving the transformation from the (X , Y , Z ) 

s s s 

system to the (X , Y , Z ) system, 
e e e 

Since the objective of the total model is to determine the position and 
orientation of a vehicle at any time in flight in different coordinate systems, 
a procedure for transforming the coordinates of a vector from one triad to 
another is mandatory. This procedure can be established by developing a 
transformation between each pair of triads. A transformation is determined 
by the right-hand rotations required to carry the first triad to the same 
orientation as the second triad. Right-hand rotations are represented by 
3x3 matrices. 

Three basic right-hand rotations occur frequently throughout this paper. 
The transformations for these rotations follow: 


m 


X 


1 0 

0 cos 9 

0 -sin 0 


0 

sin 0 
cos 0 


[ 0 ] 


y 


cos 9 
0 

sin 9 


0 -sin 9 

1 0 

0 cos 9 


16 } 


z 


cos 9 sin 9 0 

-sin 9 cos 9 0 

0 ' 0 i 


Here, [0]^ is the transformation for a right-hand rotation about the X axis 
through an angle 6 ; [0]^ is the transformation for a right-hand rotation 
about the Y axis through an angle 9 ; and [0] is the transformation for a 
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right-hand rotation about the Z axis through an angle 9 . An inspection of 

Figure 7 makes a derivation of these transformations trivial. All three 

T T 

transformation are orthogonal since [0], [0].=[0].[0] = [I] , for 

i = x, y, z where [I] is the identity matrix. Hence, the transpose equals 
the inverse. This fact will be used extensively throughout this paper. All 
of the transformations required by the model will be combinations of these 
three basic ones. 

Two rotations are required to align the surface triad with the equatorial 


triad. A right-hand rotation about the X axis through the angle A will 

s z 


carry the Z axis into the launch meridian plane (Figure 5) . Then a negative 
s 

right-hand rotation about the new Y g axis through the angle - cp^j will 
complete the alignment. Hence, 



Since the surface triad and the earth -centered triad have the same orientation, 
the above transformation will apply there also 0 This is 
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(2-3) 






and 



(2-4) 


An inspection of Figure 6 shows that the (X , Y , Z ) triad can be 

G 6 G 

aligned with the (U , V , W ) triad by these same two rotations. This 

G G G 

means that 




(2-5) 


and 



Only one rotation is required to align the (U, V, W) triad with the 
(U , V^, W^) triad (Figure 8). Thus, 





u 



u 

e 

(t> © < 

1 

= [v] 

L J x 

V 

w 

= Ltl u 

L J X 

(2-7) 


and 


u = 



( 2 - 8 ) 


The necessary transformations are now available to develop a trans- 
formation from the surface triad to the launch-point triad. From equations 
(2-1) , (2-6) , and (2-7) , 
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ROTATION OF THE EARTH 



or 





x - [TSE] x 
e s 


(2-9) 


where 


[TSE] 



The symbol TSE was chosen as an abbreviation for transformation from the 
space-fixed surface triad to the earth-fixed launch-point triad. Since [TSE] 
is a product of orthogonal matrices, [TSE] is orthogonal. Therefore, 

3T = [TSE] T F . (2-10) 

s e 


Selection of the surface triad to simulate the inertial platform system 
requires that a transformation be established from this triad to the vehicle 
system. In addition, this transformation must be a function of the simulated 
platform gimbal readings. The three right-hand rotations required to revolve 
the inertial platform triad to the same orientation as the vehicle triad were 
given in the discussion of Eulerian angles in Chapter I. Hence, the transfor- 
mation in question is a product of the transformations of those three right-hand 
rotations. That is, 


x 

m 



z 

m 


Nx M.N 


X 

y 


[TSM] x 
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where 


[TSM] - M x [, y ]J, p]y • 

Since the (X , Y , Z ) triad has the same orientation as the (X, Y, Z) 
s s s 

triad, 

x - [TSM] x ( 2-1 1) 

m s 

and 

x = [TSM] T ST . (2-12) 

s m 

This is the last transformation required by the total model. 


23 


CHAPTER HI 


EQUATIONS OF MOTION 

The purpose of the total mathematical model is to determine the 
relative motion of a Saturn vehicle during the first stage of powered flight. 

The model is required to predict the vehicle's translational motion with respect 
to the launch point, the space -fixed position and velocity which the launch 
vehicle digital computer would compute, and the three Eulerian angles which 
the vehicle's inertial platform would measure. To accomplish the first two 
objectives, all translational motion of the vehicle is first referred to the space- 
fixed surface triad. Next, translational motion of the earth-fixed launch point 
triad is determined with respect to the surface triad. Earth-fixed motion of 
the vehicle is then obtained by algebraically combining the space-fixed trans- 
lational motion of the launch-point triad with the space -fixed translation 
motion of the vehicle. To predict the three Eulerian angles, the angular 
motion of the vehicle is first computed with respect to the vehicle -fixed triad. 
This motion is transformed to Eulerian angular motion and integrated to 
obtain Eulerian angles. 

Since most of the equations of motion involve vectorial quantities such 
as position, velocity, and acceleration, vectorial notation is used extensively 
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in this chapter. A subscript is attached to each vector to indicate the triad 


to which the vector components refer. The subscripts s, c, u, ue, e, and m 
are used to refer to the surface, earth-centered, equatorial, earth-equatorial, 
launch-point and vehicle triads, respectively. Since the surface triad and the 
earth-centered triad have the same orientation, a vector referred to either 
triad has the same components. 

The equations that describe the translational motion of the launch-point 
triad with respect to the surface triad will be derived first. Figure 9 shows 
the relative position of the two triads for arbitrary t (the time from inertial 
platform release) . Notice that 


R(t)| = iR(t-O) | - R (ip ) 

o 


since the launch point remains on the Fischer Ellipsoid at a fixed geocentric 
latitude ip^ . The following three equations are obtained from an inspection 
of Figure 9: 


R <V 


cos 2 ip 

o 



1 

2 


(3-1) 


RP = R(t) - R(t=0) 
s s s 


R(t) 


ue 


R( ip ) sin ip 
o o 


-R( ip ) cos ip 
o Y o 


- 1 ue 


(3-2) 


(3-3) 
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By equations (2-2), (2-8), and (3-3) 


E(t) , = [ A JJ [% - i] y T [“„*] 


T 

'X 


R(^ ) sin ip 
o o 


-R( ip ) cos ip 
o o 


ue 


- kj; k - f] 


T 

y 


R( ip ) sin ip 
o o 

R( ip ) cos ip sin(cu t) 
o o e 

“R( ip ) cos ip cos (to t)| 

-^ue 


Then, substituting equation (3-4) into equation (3-2) yields 


(3-4) 


Rp s - N X T h - f] 


T 

y 



sin ip 


0 

R( V 

cos ip sin(co t) 
T o e 


-cos ip cos(co t) 
o e 


ue 


E <V 


sin ip 


o 


0 

y 

i 

-cos ip 


_ o_ 

ue ' 


■ n; k - f 


1 T 
J y 


R(^ ) cos ip sin(to t) 
o o e 

R (ip ) cos ip (1 - cos co t) 
'o' o e 


ue 


( 3-5) 
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T T 

Since the product A I \ C P 0 “ “|~j i s constant with respect to time, 

L J x L Jy 


RP s =■ i (RP s> “ [\T|Vt] 


c o R(ip ) cos ip cos ( co t) 
e o o e 


oo R(ip ) cos ip sin(co t) 
e o o e 


^ue 


(3-6) 

Equations (3-5) and (3-6) above describe the translational motion of the 
launch-point triad with respect to the surface triad, as required. 

The first law of motion of Newtonian mechanics is used to describe the 
translational motion of the vehicle center of gravity. This law states that a 
particle of mass m , subject to a force F , moves relative to a basic frame 
of reference in accordance with the equation 


F = km A 


(3-7) 


where A is the acceleration of the particle and k is the universal positive 
constant that depends on the choice of units of force, mass, length, and time. 
In the model being described, it is assumed that a system of units is used for 
which k = 1 . The frame of reference to which motion of the vehicle CG is 
referred is the space-fixed surface triad. 

Notice that neither m nor F is constant in the model. Mass is a 
function of the rate at which the vehicle's engines consume the stored propel- 
lant. F depends on the magnitude and direction of the engines' thrust, local 
atmospheric conditions, velocity of the vehicle with respect to the air, and 
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attitude of the vehicle. The forces that act upon the vehicle are described in 
Chapter IV. 

Equation (3-7) gives, on resolution into components along the surface 
triad axes, 


m 


d 2 X 

S 

dt 2 


F 

xs 


d 2 Y 


m 


dt 2 


F 

ys 


d 2 Z 


m 


dr 


F 

zs 


where F , F , and F are the components of force along the axes, 
xs ys zs 

These three equations are each solved for acceleration and integrated to obtain 
velocity. That is, for arbitrary t , 


dX 
s 

dt 

dY 
s_ 

dt 

and 

dZ 
s 

dt 



(3-8) 


(3-9) 


(3-10) 


Position of the vehicle at arbitrary t is determined by integrating 
equations (3-8), (3-9), and (3-10). c 



(3-11) 
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dt 


(3-12) 



(3-13) 


For a length of time A after the inertial platform is released, a 
Saturn vehicle does not move from its launch pad. For this reason, a second 
time reference, r = t - A , is defined. If r < 0 , the vehicle CG has a 
constant velocity because of the earth’s rotation. Figure 10 shows the relative 
position of the vehicle for 0 < t < A . During this time, 


RS 


RP 


[tse] 


REO 


(3-14) 


and 


RS = -37 ( RS ) = ( co ) x RC 
s dt s V eJ s 

s 


(3-15) 


where 


— v r 

T 

. i t r 7 t 1 

00 

e 

(“a) = L 

4 z k - t 

0 

x s 

L 

X 

f 

1 

'C 

0 


(3-16) 
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and 


RC = R(t) + [TSE] T REO . (3-17) 

s s e 

Equations (3-14) and (3-15) are not valid for r > 0 . Figure 11 shows 
the relative position of the CG after first motion. At this time, 


RE = [TSE] (RS - RP ) 
e s s 


( 3-18) 


and 


RS 


s 



+ [TSE] T RE e 


(3-19) 


By differentiating equation (3-19) , an expression is obtained that contains the 
velocity of the CG with respect to the launch point triad. 






d T ■ 

,* |tse! i 


RE + [TSE] 


4 R E 

dt e 


d -~— 

This equation can be solved for — ( RE^) by transposing terms and multi- 
plying each one by [TSE] . 



[TSE] ^(RS g ~ RP g ) - [TSE] (^[TSE] T jRE e 
[TSE](RS g - RP g ) - [TSE] ^ [TSE] T ') RE g 


(3-20) 


All the terms in equation ( 3-20) have been previously defined except 
d T 

— [TSE] . An expression for this term will now be derived, 
dt 
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Define I , J , K to be unit vectors along the X , 
s s s s 

I , J , K to be unit vectors along the X , Y , Z axes, 
e e e e e e 

(2-10) , these two unit triads are related by the equation 


Y , Z axes and 
s s 

By equation 


I 


I 

s 


e 

J 

= [TSE] T 

J 

s 

e 

K 


K 

s 


e 


Notice that the derivatives of the unit vectors I , J , K are nonzero since 

e e e 

their directions continuously change because of the earth ! s rotational velocity. 
Hence, 


dr? 

dt 



V 


(3-22) 


for 77 — 1 , J , K where 
' e e e 



CO 


el 


co 


e2 


co 


e3 



(3-23) 


By defining [TSE] = [A ] where i, j = 1,2,3, equation (3-21) can be 
rewritten in component form as 


! g = A U I e + A 21 J e + A 31 K e 
J g = A 12 Iq + A 22 J e + A 32 K e 

K g = A 13 I e + A 23 J e + A 33 K e 
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An application of equation (3-22) and the fact that the unit vectors I , J , and 

s s 

K do not vary with time yield expressions that contain the derivatives of all 
s 

the elements of [TSE] . 



(3-25) 


and 

dK 

dt = 0 = A 13 I + A 23 J g + A 33 K e + x (A 13 I + A 23 J e + A 33 K e ) , 

(3-26) 

where 

d . 

A.. = — A.. 

ij dt i] 

By carrying out the indicated cross product and collecting terms, equation 
(3-24) becomes 

0 = (An + w g2 A 31 - aj e3 A 2i) + ( A 2i + Cc) e 3 A ii “ “el Aai ^ J e 

+ (A S1 + o) ei A 21 - o> e 2 A 11 ) K g . (3-27) 

The vectors I , J , and K are linearly independent since they are mutually 
e e e 

orthogonal. Hence, 
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and 


A 11 = co e 3 A 2i " w e 2 A 3i » 


A 21 = A 31 - cu _ 0 A 


1 "1 ~ 


A 3i = u e 2 A n ~ W ei A21 1 ' 

An application of the above procedure to equations (3-25) and (3-26) yields 
expressions for A., where i = 1,2,3 and j = 2 , 3 . These expressions are: 


A 12 = 

= CD 

e3 

A 22 

w e2 

A 32 

A 22 = 

= w el 

A 32 

" “e3 

A 12 

A 32 = 

W e2 

A 12 

~ w ei 

A 22 

A 13 = 

= "e3 

A 23 

W e2 

A 33 

A 23 = 

= w el 

A 33 

- CD « 
e3 

A 13 

A33 = 

W e2 

A 13 

' CD 

el 

A 23 


Using the above equations, we can now express the time derivative of 
[TSE] as 


dt 


[TSE] = 


11 

A 12 

A 13 

21 

A 22 

A 23 

31 

A 32 

A 33 
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- " el A 2 i 

CO 
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~ U e3 
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co A 
el 


A 2 i 

A 22 

A 23 



W e2 

-co , 
el 
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A 31 

A 32 

A 33 




where 



= [*] 

[TSE] 

3 


0 

CO 

e3 

~ 0; e2 

*] = 

" W e3 

0 

“el 


W e2 

-co 

el 

0 


(3-28) 


Therefore, 


dt 


[TSE] 


(^ ITSE] ) = ([^][TSE]) T = [TSE] T [$] T 


(3-29) 


By substituting equation (3-29) into equation (3-20), — (RE ) - RE can be 

d t e e 

be expressed as 


RE - [ TSE] ( RS -RP ) - [*] RE 

e s s e 


= [TSE] , RS s - RP s ) - (m 6 ) 

x ' e 


x RE 


(3-30) 
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Equations (3-18) and (3-30) describe the relative translational motion of the 
vehicle CG with respect to the launch-point triad as required. 

Total motion of a vehicle is described by the model in terms of the 
translational motion of its CG and rotational motion of the vehicle about its 
CG. The equations that describe the second type of motion are all that remain 
to be derived. 

The model assumes a vehicle to be a rigid body constrained to rotate 
about its CG. Thus, the angular momentum about the CG is 

L = Icoi + Icoj + Icok , (3-31) 

x x y y z z 

where i, i, k are unit vectors along the X , Y , Z axes, respectively; 

J & m m m 

1,1,1 are the principal moments of inertia about these axes; and 
x y z 

co , co , co are the components of angular velocity co of the vehicle about 
x y 5 z 

these axes. Here, To is the angular velocity of the vehicle triad with respect 
to the surface triad. Equation (3-31) is classical and can be found in most 
books on Newtonian mechanics. 

The total moment of the external forces about the CG is 
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where 


dl 

I = — 2- 
r? dt ’ 

dw 

V 

co = ~n — 

7] dt 

for 77 = x,y, z 

Since the time rate of change of the principal moments of inertia is very 
small, it is assumed that i = 0 for rj = x,y, z . Each of the unit vectors 
i, j,k changes only in direction so that 

= "m * 5 for ' ' *' i,k ' 

Thus, equation (3-32) can now be simplified to 

M = Icoi + Icoi+Icok + co x (Icoi + Icoj+Icok) 
xx yy zz e xx yy zz 

- Tl CO - CO CO (I - I )1 i + r I CO - co co (I “I )1 j 

[xx y Z y z'J Lyy xzz xJ J 

+ I CO - CO CO (I —I ) I k 
L z z x y x yj 

This equation can be rewritten in component form as 


M t 

- I 

CO 

- CO 

CO (I 

-I ) 


X 

X 

y 

z y 

z 

m 2 

= I 

CO 

- CO 

CO (I 

-I ) 


y 

y 

X 

z z 

y 

m 3 

= i 

CO 

- CO 

CO (I 

-i ) 


z 

z 

X 

y x 

y 
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where M 1# M 2 , M 3 are the components of M along i,j,k . Therefore, 

(3-33) 
m 

An expression for the total moment vector, M , about the vehicle CG is given 
in Chapter IV. Angular velocity of a vehicle about its CG is determined at 
arbitrary t by integrating equation (3-33) to obtain 


Mj/l -h co co (I -I )/l 
x y z y z x 

M 2 /I + co (o (I -I )/l 
y x z' z y y 

M/I + co Ct> (I - I )/l 
z x y X y z 



Equations ( 3-33) and ( 3-34) describe the rotational motion of a vehicle about 
its CG as required. 

Since the orientation of a Saturn vehicle is defined by three Eulerian 
angles, some method of computing these angles must be established. One of 
the simplest methods is to transform the vehicle angular velocity to Eulerian 
angular velocity and integrate the Eulerian angular velocity to obtain Eulerian 
angles. To this end, a transformation from the vehicle angular velocity to 
Eulerian angular velocity will now be derived. 

The required transformation will be derived from infinitesimal rota- 
tions of a vehicle about its CG . Suppose that a rigid vehicle has turned through 
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an infinitesimal angle about an axis through its CG. This infinitesimal rota- 
tion may be represented by an infinitesimal vector 6 77 (Figure 12) . 6 77 

determines the displacements of all points of the vehicle. 

Let r be the position vector relative to the CG of a point A before 
the displacement and r + 6r the position vector after displacement. In the 
rotation <5 77 the point A moves perpendicular to the plane containing the 
vectors r and 6 77 , in the same direction as the vector product 677 x r . 

The magnitude of dr is \ 6 r)\AB= I <5 77 | | r | sin 0 (Figure 12). Hence, 

6r = 6 77 x r „ (3-35) 

Infinitesimal rotations about a point can be added vectorially. Let 677 
and 6m be two infinitesimal rotations applied in succession about axes 
through the vehicle CG. Then the position vector of a point A in the vehicle 
is r before displacement, r + 677X r after the 677 rotation, and 
r + 677 x r + 6m x (r + 6 77 x r) after the 6m rotation. By neglecting 
infinitesimals of the second order, the resultant displacement of point A is 

677 x r + 6m x r = (677 + 6m) x r * (3-36) 

Thus, an infinitesimal rotation of a vehicle about its CG may be described 

either as an infinitesimal rotation 677 or by means of infinitesimal increments 

in the Euleriaft angles . The rotation of a vehicle from the orientation 

(cp , <p , cp ) to a second orientation ( cp + A cp , cp +A cp , <p + Acp ) 
p y r P p y y r ^r 

may be accomplished by applying the following finite rotations in order: 
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CG r 

FIG. 12. INFINITESIMAL 


ROTATION 
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Acp^i , A K r , A cp^ J where i, K', and J are unit vectors along the 
, Z ! , and Y axes, respectively. (See Figure 3 of Chapter I. ) If the 
increments in the Eulerian angles are infinitesimal, the order of these rota- 
tions is immaterial, so that 

<5 Tj — A cp ^ i + A cp J 4- A cp K* . ( 3— 37) 

The infinitesimal rotation of concern here is the rotation 173 I dt in the 

direction of the unit vector 73/1731 where 73 is the angular velocity vector 

of the vehicle. Corresponding to this rotation, ( 1 73 I dt) (73 / 173 I ) = 73 dt , 

are the three infinitesimal increments to the Eulerian angles, dcp i, dtp J, 

r p 

dcPyK 1 . Thus, equation (3-37) becomes 

co dt — dc^i + d^J + dt^K. 1 ■ (3—38) 

By dividing each term of equation (3-38) by dt , the expression 

co - i i + i J + i K' (3-39) 

r ^p y 

d cp 

. ' m 

is obtained where cp = — for m = r, p,y . It is now convenient to 

^m dt J 

convert equation ( 3-39) to an expression involving unit vectors along the 
X , Y , and Z axes of Figure 3. Define j , k , I f , J T , I , and K to be unit 
vectors along the , Z^ , X T , Y’ , X, and Z axes, respectively. Then, 
by an inspection of Figure 3, we see 
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(3-40) 


I' 

J 

‘ N y 

I 

J 

_K’_ 


_K_ 


and 

I'" 

J 

K'_ 

From these two equations, it is easy to show that 

i m cos cp cos cp I + sin cp J - cos cp sin cp K 



and 


(3-41) 


(3-42) 


K T = sin <p I + cos cp K t (3-43) 

P P 

Then, by substituting equations (3-42) and (3-43) into equation (3-39) and 
simplifying, we obtain 

To ~ (cp cos cp cos cp + <p sin cp )I + ( cp sin cp + cp ) J 
s r r y y r p ^r ^y v p 

+ (<p cos cp - cp cos cp sin cp ) K . (3-44) 

y P r y P 

A resolution of a; in the vehicle system can now be obtained by using 
equation ( 2-11) : 
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oo = [TSM] 
m 


' cp cos cp cos cp + cp sin cp 
r ^p 

cp sin cp + <b 

r Y y 

cp cos cp - cp cos cp sin cp 

l y p r y p. 


where 


cp sin cp + cp 
p y r 

cp sin cp + cp cos cp cos cp 
r y r ' p ^r ^y 

cp cos cp - <p sin cp cos cp 
r y r p r ^y 


sin cp 

P 

0 

1 


1 

a 

•S- 

L 

cos cp cos cp 
r r y 

sin cp 

r 

0 


<p 

y 

-sin cp cos cp 
L r y 

cos cp 

r 

0 


<p 

L r J 


[TEM] 


?> 




</> 


( 3 - 45 ) 


[TEM] = 


sin cp 0 


cos cp cos cp sin cp 
r y i 


-sin cp cos cp cos cp 
r y i 


i 

0 

0 


Here, the symbol TEM was chosen to stand for the transformation from 
Eulerian angular rates to vehicle angular rates. The inverse of [TEM] can 
be computed by multiplying the transpose of the matrix of cofactors of [TEM] 
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by the scalar 1/ I [TEM] I , where I [TEM] 1 is the value of the determinant 
of [TEM] . Thus, 


[TEM] 


-1 


1 


cos cp 


0 

cos cp 

r 

-sin cp 


sin cp cos cp 

r y 


cos cp cos cp 

r y 


cos cp 


-sin cp cos cp 
y r 

sin cp sin cp 

y r _ 


0 cos cp sec cp 
r y 

0 sin cp 


1 -tan cp cos cp 

. y r 


-sin w sec cp 
r y 

COS <£> 


tan sin <vp 
y r 


(3-46) 


Therefore, [TEM] 1 exists when cp ± ±90 degrees, so that 


cp 

p 


CO 

X 




= [TEM] -1 

CO 

y 

= [TEM] -1 73 

m 

(3-47) 

A- 


CO 

L z J 




Since the yaw gimbal of a Saturn inertial platform is physically restricted to 
rotations of ±45 degrees from its initial position, the above restrictions on 
cp create no problem. 

y 

Equation ( 3-47) is the required transformation for converting the 
vehicle angular velocity to Eulerian angular velocity. The Eulerian angles 
cp^ 9 ^y anC * ^r are com P u ^ e( ^ arbitrary t as follows: 


<P 


V 




(3-48) 
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for tj — p, y, r . This is the last equation required to describe the rotational 
motion of a rigid vehicle about its CG. 



CHAPTER IV 


FORCE AND MOMENT EQUATIONS 

A Saturn vehicle experiences forces from three basic sources during 
its first stage of powered flight. These sources are the vehicles first-stage 
engines and the earth T s gravitational potential and atmosphere. This chapter 
provides the equations required to compute the total force and moment vectors 
that result from these forces. 

The S-IC stage of a Saturn V vehicle has five engines (Figure 1) that 
propel the vehicle through the earth T s atmosphere. The four outboard engines 
are attached to the base of the stage by mechanical joints that permit the 
engines to gimbal or rotate through a limited angle from their null position. 
Gimbaling of these four engines according to commands from the flight control 
computer forces the vehicle to follow the desired attitude commands. Each 
outboard engine has two hydraulic actuators, attached 90 degrees apart, that 
execute the flight control computer commands. One actuator provides 
attitude control in the pitch plane, and the other actuator provides attitude 
control in the yaw plane. Attitude control is achieved in the roll plane by 
moving both the pitch and yaw actuators in the proper directions. 
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In Chapter I, the word "control” was defined, and the control system 
was briefly described. A more detailed description of this system will accrue 
as the model for it is presented. 

A block diagram of the control system for a Saturn vehicle is shown in 
Figure 13. Attitude error signals are computed in the LVDC from the plat- 
form gimbal angles and the commanded attitude angles. The vehicle angular 

velocities about the X , Y , and Z axes are measured by three rate 

m m m 

gyroscopes. In the flight control computer, signals from the rate gyroscopes 
and the attitude error signals are electrically filtered, multiplied by a gain 
factor, and then combined to generate the commands to the four outboard 
engine actuators. 

The commanded attitude angles x > X^ , and X r are defined to be the 

three Eulerian angles through which the inertial platform triad must be 

rotated to align it with the desired orientation of the vehicle. If i , , and 

k are unit vectors along the vehicle X , Y , and Z axes, respectively, 
c m m m 

when the vehicle has the orientation ( x^ , Xy > X r ) > then 
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j c 
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L c - 

■ M x W z [v] y 

1 ft c-t 

m ® 

1 

= [TSC] 

J 

s 
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s__ 


where [TSC] - ["x 1 |~X 1 [x 1 and I , J g , and K g are unit vectors 

L r J x L *- ^y 

along the surface triad X , Y , and Z axes respectively. 
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Similarly, the attitude error signals ^ , ip ^ , and ip^ are defined to be the 

three Eulerian angles through which the vehicle triad (X^ # , ZJ must 

be rotated to align it with the (i , j , k ) triad. Therefore, 

c c c 




i 

[VLM.M, 

j 


k 


(4-2) 


where i, i, and k are unit vectors along the vehicle X , Y , and Z 

J m m m 

axes respectively. By equations (2-12) and (4-1) , a second expression for the 

transformation from the (X , Y , Z ) triad to the (i , j , k ) triad is 

m m m c c c 

obtained. 


" i 
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j c 

= [TSC]ITSM] T 

i 

3 

k 
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Hence, by equations (4-2) and (4-3), 


p r ] [tf y ] [* P ] = [TSC][TSM] T . (4-4) 

The following expressions for ip ip , and ip^ can be obtained from 
equation (4-4) by performing the indicated product and equating corresponding 


elements of the two matrices. 



ip = (<p - x ) cos 
y y y 


+ x ' 

r r 


- ( <P p - Xp) cos 


+ X 

_y y 


sm 


cp + X 
r r 




(4-6) 


^ = (<p-X)+(^“X)sin 

r r r p p 


(p + X ' 

ly x 


(4-7) 


The details involved in deriving the above expressions and the approximations 
on which they are. based are given in Reference 3. These are the expressions 
used by the LVDC for computing the attitude errors. 

The attitude errors are converted to analogue signals by the LVDA 
before they reach the flight control computer. Here, both the attitude error 
signals and the angular velocity signals from the rate gyroscopes are passed 
through electrical networks called filters. These networks are designed to 
remove the effects of flexible body motion from the input signals and to main- 
tain proper control system stability. In this paper the symbols ip^ 9 ip , ip^ 
and 0 ^ , are used to denote the filtered attitude error and angular 

velocity signals, respectively. 

A special digital filtering technique (Reference 4) has been developed 
to simulate the electrical control filters. This technique will not be discussed 
in this paper since it is widely used and a discussion of it is lengthy. 

The filtered attitude error and angular velocity signals are combined 
in the flight control computer by the following control laws to give the 
commanded thrust deflections in the pitch, yaw, and roll planes. 
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Pitch: 


(4-8) 


Yaw: 


Roll: 


P = A ip + A <p 
pc Op p lp ^p 


P = A n ip + A a <f> 
yc Oy y ly y 


B = A ip + A. (h 
rc Or r lr 


(4-9) 

(4-10) 


Here, and for j] = p,y, r are gain factors that provide stability to 

the control system. The magnitudes of these factors are specified as step 
functions of flight time. 

Individual commands to the actuators of the four control engines are 
defined by the following equations: 
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(4-14) 

(4-15) 

(4-16) 

(4-17) 
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where p and p are the commands to the pitch and yaw actuators of the 
i th outboard engine. 

The positive signs for all the control parameters are shown by Figure 

1 . A positive attitude error about a vehicle axis requires a negative right- 

hand rotation about that axis to correct for it. A positive angular velocity 

about a vehicle axis has the same sense as a positive right-hand rotation about 

that axis. The bases of the four outboard engines are moved parallel to the 

+Z ffi axis by a +p ^ creating a negative moment about the pitch or axis. 

Similarly, a +p ^ moves the base of the engines parallel to the negative 

axis, creating a negative moment about the yaw or axis. Thus, positive 

actuator movements correct for positive attitude errors. 

The thrust vector of each of the five S-IC stage engines is ideally 

parallel to the axis when the engine is in the null position. However, an 

engine’s thrust vector alignment may deviate because of electrical and 

mechanical tolerances. To simulate the effect of the thrust vector deviation, 

the null position of each engine is assumed to be offset by the two small angles 

A d . and Ad ., for i = i,2,3, 4, 5. Here, A d .and Ad . have the same 
K pi yi pi yi 

sense as B and p in Figure 1. 

p y 

The method used to resolve the components of force of an S-IC stage 

engine along the vehicle axis is shown in Figure 14. In this figure, the 

(X M , Y" , Z TT ) triad has the same orientation as the (X , Y , Z ) triad. 

m m m 

F^. is the thrust vector of engine i for i = 1, 2, 3, 4, 5 . Note that 

ft = p = 0 since the center engine does not gimbal (Figure 1) . Let F . , 
po yo exi 
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F . . and F . be the components of F . along the X" , Y” , Z n axes, 
eyi ezi ei 

respectively. From an inspection of Figure 14, it is easy to show that 

F . = IF . I/K (4-19) 

exi ei 

F . = |F .1 tan (j8 . + A/J .)/K (4-20) 

eyi ei r yi ^ yi 

and 

F . = - |F . I tan (P . + A/3 .)/K (4-21) 

ezi ei pi pi 

where 

K = n/ 1 + tan 2 (jS ~ + A/3 .) + tan * ((3 . + A/3 .) 

pi pi yi yi 

The thrust of each S-IC stage engine is a function of time and atmos- 
pheric pressure at the engine. |F .1 is computed from the expression 

IF . I = F . - P A A. , (4-22) 

ei vi A i 

where F . is the thrust of engine i firing in a vacuum, P. is the atmos- 
vi A 

pheric pressure at the engine, and A is the exit area of engine i . F^. is 

predicted as a function of time. 

As the vehicles engines consume propellant, the vehicle CG and thus 
the origin of the (X . Y , Z ) triad change with respect to the vehicle 
frame. A frame-fixed triad (X , Y f , Z f ) is defined to serve as a reference 
for referring the location of the CG, the attach points of the engines, and 
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certain aerodynamic characteristics. The origin of this triad is at the inter- 
section of the geometrical centerline of the vehicle and the plane that contains 
the attach points of the S-IC stage engines. The +X f axis is directed along 
the stage centerline, and the Y f and Z f axes are parallel to the and 
axes, respectively. Thus, this is a right-hand triad with the same orientation 

as the (X , Y , Z ) triad (Figure 15). 
m m m 

An expression for the moment of an S-IC engine about the vehicle CG 
can now be derived. The total thrust of an engine is assumed to be applied 
at the attach point of the engine to the stage. Figure 15 shows the relative 
position of the vehicle CG with respect to the attach point of engine i . By 
definition, the moment M of a vector V about a point 0 is 

M = ? x V , (4-23) 

where r is the vector directed from point 0 to the point of application of V . 
Therefore, from Figure 15, the moment JML of the thrust of engine i about 
the vehicle CG is 

M. = r. x F . = (P. - CG) x F . , (4-24) 

i i ei i ei 

where and CG are the frame-fixed position vectors of the attach point of 
engine i and the vehicle CG, respectively. 

The above equation along with equations (4-19) , (4-20) , and (4-21) 
provide the total model with expressions for computing the forces and moments 
resulting from the thrust of the S-IC stage engines. 
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A Newtonian potential function (Reference 1) is used to determine the 
earth’s gravitational acceleration vector at the vehicle CG. This function was 
discussed briefly in Chapter I. No further discussion of this function will be 
presented here except to say that it uses the position vector, PP^, °f the 
vehicle CG referred to the (U , V , W) equatorial triad to compute the 
gravitational acceleration of the vehicle. 

The position vector PP^ will be computed from the vector PP c 
(position of the vehicle CG relative to the earth-centered triad) . From 
Figures 9 and 11, 

PP - R(t=0) + RS . (4-25) 

c s s 

Then by equation ( 2-3) 



PP is used as the independent variable in the Newtonian potential function to 
obtain the gravitational acceleration of the vehicle, GR^ • 

Since the equations of motion refer all translational motion of a vehicle 
to the surface triad, it is convenient to refer the gravitational acceleration to 
it also. Thus, by equation (2-2) 


GR 

s 



(4-27) 
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The gravitational force vector, FG , referred to the surface triad, can now 

s 

be derived from an application of equation (3-7) , the first law of motion of 
Newtonian Mechanics. Hence, 

FG = mGR * (4-28) 

s s 

where m is the mass of the vehicle. 

As a vehicle flies through the earth’s atmosphere, it experiences the 
forces of buoyance and air friction. The buoyance force, F^ , is assumed to 
have the same direction as the position vector of the vehicle relative to the 
center of the earth. Its magnitude is 

|F I = V p |GR I , (4-29) 

b A s 

where V is the volume of air displaced by the vehicle and p is the density 
of the air through which the vehicle is flying. Therefore, 

(\\ ’ V A» |5 V PV l?5 c' ■ (4 - 30) 

Notice that PP / |PP I is a unit vector that has the same direction as F . . 

c c b 

The point of application of F^ is assumed to be at the vehicle T s CG. 

Therefore, the moment of F, about the CG is zero. 

b 

Figure 16 shows a vehicle. moving through the earth’s atmosphere with 
velocity RV with respect to the surrounding air. Contact of the vehicle with 
particles of air produces a distribution of pressure over the vehicle’s entire 
surface. Rather than consider the pressure distribution, its effect on the 
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motion of the vehicle is accounted for by applying a suitable single force at 

a point on the vehicle’s geometric centerline called the center of pressure 

(CP) . To preserve the laws of Newtonian mechanics, the single force, F^ , 

must be equal to the vector sum of the forces of the pressure distribution, and 

the CP must be located at a point that causes the moment of F . about the 

A 

vehicle CG to equal the moment of the pressure distribution about the CG. 

Since the pressure distribution directly opposes the vehicle’s relative velocity, 

F. has the same sense as that of -RV . 

A 

The total force, F^ , is customarily resolved into two components. 

The normal force, F^ , is the component normal to the vehicle’s longitudinal 

axis and the drag, D , is the component parallel to the longitudinal axis. 

Thus, to describe the aerodynamic force of the pressure distribution, the 

variation of F , D, CP, a, ot A and RV must be specified, 
n i 

For a given Mach number and a (angle -of -attack) , F and D are 

n 

proportional to the aerodynamic pressure, Q , and the cross-sectional area 
of the missile, A . Here, 

Mach - i RV l/VS , (4-31) 

where VS is the velocity of sound in the air surrounding the vehicle and 

Q = p IRV l 2 /2 . (4-32) 

Because of these proportionalities, it is customary to write 

F - C QA (4-33) 

n n v ' 
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and 


(4-34) 


D = C , Q A 
d 

where and are dimensionless proportionality coefficients. 

Both coefficients and the distance from the origin of the (X f , Y f , Z^) 

triad to the CP are determined experimentally as functions of Mach and a . 

Resolution of F^ into drag and normal force and the moment of about 

the CG are directly dependent upon a . Dependence on Mach accounts for 

variation of these three parameters with air temperature as well as relative 

velocity. The main variation in the velocity of sound arises from the fact 

that the velocity is proportional to the square root of the absolute temperature. 

RV is the difference between the velocity of the vehicle and the velocity 

of the air surrounding the vehicle. Since the earth T s atmosphere rotates with 

the earth, the velocity of air has a component because of the rotation of the 

earth, ( oj J x PP c , and a second component because of the velocity of the 
' 's 

air with respect to the earth, W . Thus, RV referred to the (X , Y , Z ) 
^ e m’ m m 

triad is 


RV 


m 


[TSMjj RS c 


( a>e) xPP c - [TSE) T W g 


(4-35) 


The parameters a and a : (Figure 16) are determined from the 

components of RV^ • The expressions used for computing them are 


RV 


a - arc cos 


(4-36) 


IRV 


m 
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and 


o' = arc tan 


RV 
2 

RV 


( 4-37) 


where RV , RV and RV are the components of RV along the X , Y , 
x * y z m m m 

Z m axes ’ respectively. 

Expressions for the components of F^ along the vehicle-fixed axes 
can now be established. From Figure 16, it is easy to see that 



-D 


-F sin cv 
n i 


-F cos a, 
n i J 

m 


(4-38) 


where D and F^ are defined by equations ( 4-33) and ( 4-34) . 

Equation (4-23) is used to determine an expression for the moment of 

F. about the vehicle CG. From Figures 15 and 16, it is apparent that 
A 

- CG (4-39) 

m ' 

m 


CP 

0 

0 


where r* corresponds to Y in equation (4-23) and CP is the distance from 

the origin of the (X^, , Z^) triad to the CP . Therefore, the aerodynamic 

moment, M. , of F about the CG is 
A A 

(“a) X ( ^ A ) • (4 - 40 » 

\ / m \ / m 
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Since the atmospheric parameters p , air density, P^ , atmospheric 
pressure, and VS , the velocity of sound, are functions of the altitude of the 
vehicle above the earth, an expression for altitude must be derived. Figure 
17 shows a meridian through the polar axis of the earth that contains the posi- 
tion vector, PP , of the vehicle CG with respect to the (X , Y , Z ) 

’ c c c c 

triad. The following expressions are obvious from an inspection of the figure. 



( 4-41) 



h' - |PP I - c (4-43) 

c 

The variable h’ in equation (4-43) is assumed to be a close approximation 
to the actual altitude h . 

Expressions have now been developed for all forces and moments that 
a vehicle experiences as it moves through the earth’s atmosphere. Thus, 
expressions for the total force and moment vectors can now be obtained by 
algebraically combining the individual contributors. 

The total force is composed of the thrust of the five S-IC stage engines, 
the force of gravity, the buoyance force, and the aerodynamic force. Thus, 
the equation for the total force F g is 
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'-INTERSECTION 
OF MERIDIAN WITH 
EQUATORIAL PLANE 



h * h' 

|pp c | = c + h' 

EQUATORIAL RADIUS OF EARTH 
POLAR RADIUS OF EARTH 


17 . 


ALTITUDE OF VEHICLE 




F 

s 


[TSM] T 



+ FG s + ( r b) + [ ts “) T ( f a) •' 
s m 


( 4 - 44 ) 


The total moment about the vehicle CG is composed of the moments 
resulting from the five S-IC stage engines and the aerodynamic force. 
Therefore, the expression for the total moment vector, , is 


CHAPTER V 


NUMERICAL INTEGRATION 

The design of a complex and expensive vehicle system such as the 
Saturn V requires an analysis of the response of a vehicle to different designs, 
vehicle engineering tolerances, and various atmospheric conditions. Because 
a large part of the analysis does not involve flexible motion of a vehicle, the 
rigid-body response can be predicted by electronic computer programs based 
on mathematical models like the one developed in this paper. Programs 
of this type require an efficient numerical integration technique that will 
provide a solution within the accuracy of the data which describe the vehicle 
dynamic response characteristics. One such program, the liftoff program, 
was developed from the mathematical model presented in Chapters II, III, and 
IV. 

In a search for the most acceptable integration technique, the fourth- 
order Runge-Kutta formula and two of Dr. E. B. Shanks' integration formulas 
were each used in the liftoff program to compute the first stage trajectory of 
a typical Saturn V vehicle. Comparative data were generated for each formula 
over a range of integration step-sizes. Also, comparative data were 
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established for a technique associated with one of Dr. Shanks' formulas 

that regulates the integration step-size. An analysis of these data and a 

description of the integration formulas are presented in this chapter. 

The three numerical integration techniques to be compared are the 

fourth-order Runge-Kutta formula, the fourth-order formula developed on 

page 8 of Reference 5, which shall be referred to as formula 4-3, and the 

fifth -order formula 5-5 presented in Reference 6. Only a description of the 

dy 

equations as they apply to a differential equation of the form y = = f(t,y) 

need be given since the system of first and second order differential equations 
presented in Chapter III are reducible to an application of this result. 

The Runge-Kutta formula is presented first. If the initial values of 
the differential equation y T = f(t, y) are t 0 , y 0 , the value of y at t 0 + h 
is computed from the formulas 


*0 - f(to> yo) 


'(* 


ft . h hfo \ 

+ 2 ’ + 2 ) 


u - 


Ah h fj \ 

yo + 2 ’ y ° + 2 ) 

-- f(t 0 + h , y 0 + hf 2 ) 


( 5-1) 


y = y° + | (f 0 + 2fj + 2f 2 + f 3 ) 

evaluated in the given order. A derivation of this technique is given in most 
books on numerical analysis. 
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The formula 5-5 expressions for determining the value of y at 


t 0 + h are 

fo = f (t 0 > Yo) 

f t = f(t 0 + h/9000 , y 0 + hf 0 /9OOO) 

f 2 = f(t 0 + 3h/ 1 0. , y 0 - 404. 7hf 0 + 405hf 1 ) 

f 3 = f(t 0 + 3h/4 , y 0 + 2O241hf 0 /8 - 20250hf 1 /8 + 15hf z /8) (5-2) 

f 4 = f(t 0 + h , y 0 - 931O41hf 0 /81 + 93150011^/81 
- 490hf 2 /81 + 112hf a /8l) 

Y = Yo + ( 1 05f 0 + 500f 2 + 448f 3 + 81f 4 ) 

where each formula is evaluated in the given order. 

Formula 4-3 differs from the above formulas in that it uses 
f[t 0 - h , y(t 0 - h) ] , the value of f at t 0 - h,in computing y. The value of 
y at t 0 + h is determined from the expressions 

fo = f(t 0 , Yo) 

fi = f(t 0 -h , y(t 0 - h)J 

f 2 = f[t 0 + h/2, y 0 + h( 5f 0 - f 4 ) / 8] (5-3) 

f 3 = f[t 0 + h j Yo + h(-3f 0 + fj + 4f 2 )/2] 

Y = Yo + h(fo + ^2 + fs)/® 
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The fourth-order Runge-Kutta formula is used for the first integration 
step, then formula 4-3 is used to continue the integration process. 

A method of controlling the integration step-size to obtain the required 
accuracy is presented in Reference 7. A method is introduced for computing 
an estimate of the integration error made at each step. The error estimate, 
called a regulator, is used to determine the step-size for the next integration 
step. 

The regulator for formula 4-3 is given in Reference 5. By combining 
the T of equation (5-3) by the formula 

Yl - yo + h(-5f 0 + f 1+ 12f 2 -2f 3 )/6 , (5-4) 

a third-order solution for y at t 0 + h is obtained. The regulator, R , is 
defined to be the absolute value of the difference between the values of y 
computed from equations (5-3) and (5-4) ; that is, 

R = ly-yil - | |6f 0 - fj - 8f 2 + 3f 3 | . (5-5) 

Thus, R is a fourth-order estimate of the third-order error. Note that the 
number of additional calculations required to compute R are insignificant 
since the same evaluations of f are used to compute both y and R . 

When the integration process extends over a number of steps, the value 
of R can be monitored and used as an error indicator for controlling the step- 
size. To accomplish this, a lower limit and an upper limit L 2 (0<L 1 <L 2 ) 
are established for R . At the end of each integration step, R is compared 
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with these limits. If R < L* , the step-size is doubled; if R > L 2 , the step- 
size is halved; if Lj < R < L 2 , the step-size is not changed. When the 
step-size is changed, the fourth -order Runge-Kutta formula is used over the 
next step to establish a value for in equation ( 5-3) . The integration 
process is then continued with formula 4-3, 

The differential equations of the liftoff program reduce to a system of 
twelve first order differential equations. Thus, twelve distinct values of R 
can be calculated at the end of each integration step. Which of the twelve 
values of R or what combination of the twelve values of R should be moni- 
tored for controlling the integration step-size is an important question that 
must be answered in order to effectively use this procedure in the liftoff 
program. This question can be answered by determining which distinct values 
of R are most sensitive to the total accuracy of formula 4-3 in integrating 
the equations of motion. The data required to obtain this information can be 
described more easily after the comparative data for the three integration 
routines are presented. 

Before the comparative data were generated, a computer program was 
written and verified for each of the three integration formulas. The programs 
were written for the SDS-930 computer at Marshall Space Flight Center in 
Huntsville, Alabama. The differential equation x 2 y’ ? + 2xy’ + 2y-f-4 + x~0 
which has y - c t x + c 2 /x 2 - 2 + x 2 / 4 as a general solution was chosen to be 
used in verifying the programs. Choosing x = 1 , y = 1/4 , and y f = -25/8 
as initial values, the exact solution for x t 0 becomes y = x/8 + 15/8x 2 + x 2 /4. 
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I 


This equation was then numerically integrated over the interval x = . 1 to 
x = 5 to obtain the following data. 


Runge-Kutta 4-4 


Formula 4-3 Formula 5-5 


Step-Size 
. 10 
. 05 
.025 
. 01 
.005 
.0025 
. 001 
. 0005 


Error 

10 

X 

o 1 

95 

X 

10 2 

65 

X 

10 1 

16 

X 

10° 

10 

X 

10 _1 

62 

X 

10" 3 

16 

X 

i— 

o 

1 

rf*. 

30 

X 

10~ 5 


S tep-Size 
. 10 
.05 
. 025 
. 01 
.005 
. 0025 
. 001 
. 0005 


Error 


11 

X 

10 4 

10 

X 

10 s 

61 

X 

10 1 

53 

X 

10 _1 

32 

X 

10 _1 

60 

X 

10 -2 

46 

X 

10~ 3 

59 

X 

10 -4 


S tep-Size 
. 10 
. 05 
. 025 
. 01 
. 005 
. 0025 
.001 


Error 
. 10 x 10 3 
.46 x 10 1 
. 15 x 10° 

. 15 x 10~ 2 
.48 x 10~ 4 
.28 x 10 -5 
. 12 x 10“ 5 


Next, the trajectory parameters to be analyzed for comparing the 

different integration formulas in the liftoff program had to be selected. A 

time history of the vehicle position and orientation with respect to the 

(X g , Y g , Z^) triad basically defines the vehicle trajectory. For this 

reason, the parameters xs , ys , zs , <p , w , and cp are the obvious 

p p P P y r 

ones to be contrasted. 

Selection of data for use in generating the comparison results 
was the next task. Data which represent the dynamic characteristics of a 
typical Saturn V vehicle were chosen. Large dynamic disturbances were 
imposed on the trajectory by failing an S-IC stage outboard engine 70 seconds 
after liftoff and disturbing the atmosphere with a wind that peaked immediately 
following the engine failure (Figure 18) . The wind was directed normal to the 
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vehicle flight path in the direction that would create an aerodynamic moment 
on the vehicle having the same sense as the moment created by the engine 
failure. When the entire wind profile is increased by 1 meter-per-second, the 
control system is incapable of continuing to control the vehicle's attitude after 
the engine failure. That is, the aerodynamic moment forces the vehicle to 
begin to tumble. Thus, the chosen trajectory is representative of nominal 
flight before 70 seconds and highly disturbed flight thereafter. 

The liftoff program is written in single -precision for a CDC-3200 
computer at the Marshall Space Flight Center in Huntsville, Alabama. 

(Single -precision on this computer means that eleven digits other than the 
exponent are assigned to each floating point number. ) Tables 1, 2, and 3 
present comparative data obtained from this program for different integration 
formulas and step-sizes at the 70, 80, and 90 second time points of the 
trajectory. The three different time points in the trajectory are presented 
in case any one formula is biased at any one of the time points. 

The data presented in the tables were obtained from an on-line printer 
while the trajectories were being computed. Since the CDC-3200 computer 
does not print and compute simultaneously, run time is affected by the amount 
of output obtained. For this reason, only a small amount of output was 
obtained from the trajectories for which the computation times are presented 
in Tables 1, 2, and 3. During the first 70 seconds of each trajectory, output 
was obtained every 10 seconds; for the remainder of the trajectory output 
was obtained every 5 seconds. 
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Some way of comparing the accuracy of the solutions presented in the 
tables must be established in order to evaluate the comparative data. At 
present, no means of determining an exact solution to a trajectory based on a 
specific set of vehicle dynamic data exist. A review of Tables 1, 2, and 3 
reveals that the differences between the solutions reached by the three formu- 
las are less than 1 percent in all parameters except the roll and yaw attitude 
angles which differ less than 4 and 6 percent, respectively, at the 90 second 
time point. Therefore, it shall be assumed that the exact solution does not 
differ more than these percentages from the solution reached by any one of the 
formulas. Based on this assumption, the accuracy of each formula can be 
judged by whether the difference between the solution it reached and the 
solutions reached by the other two formulas is within the accuracy of the data 
that describe the dynamic characteristics of the vehicle. 

Because of hardware and engineering tolerances, the dynamic charac- 
teristics of a vehicle cannot be predicted exactly. Publications that present 
the dynamic data for a particular Saturn V vehicle give the statistical varia- 
tions of the data. The one -sigma data tolerances are listed in Table 4. The 
effects of these tolerances on the trajectory when applied to the dynamic data 
used to generate Tables 1, 2, and 3 were assessed using formula 4-3 with an 

integration step-size of . 125. Tables 5, 6, and 7 give the solution of the tra- 
jectory at 70, 80, and 90 seconds, respectively, with the specific tolerances 

*;pplied in both the most helpful and most harmful way. Here, the most helpful 
way is defined to be an application of the tolerance that reduces the effect of the 


76 



engine failure and wind disturbances; the most harmful way is defined to be an 
application of the tolerance that increases the effect of the engine failure and wind 
disturbances. The effect of each tolerance can be assessed by comparing the 
formula 4-3 solution of Tables i, 2, and 3 with Tables 5, 6, and 7. 

A comparison of Tables 5, 6, and 7 with Tables 1, 2, and 3 shows 
that differences in the solutions reached by all three integration techniques 
are within the accuracy of the Saturn V dynamic data. Therefore, by this 
method of judging accuracy, all three integration techniques are acceptable 
from an accuracy stand point. 

A study of Tables 1, 2, and 3 reveals that formula 4-3 would be the 
most desirable integration formula to use in the liftoff program from a stand- 
point of computer run time, accuracy, and stability of the solution. At all 
three time points, the three formulas reach a stable solution at the same 
step-size. Also, the solution reached by formula 4-3 does not differ signifi- 
cantly from the solutions reached by the other formulas, and this difference 
is within the accuracy of the input data. The differences between the Runge- 
Kutta solution and the formula 4-3 solution are so small that the solutions are, 
for all practical purposes, identical. This fact is important since the fourth- 
order Runge-Kutta formula is presently being used almost exclusively in all 
digital trajectory programs at Marshall Space Flight Center. Computer time 
would be saved by use of the formula 4-3 since it requires less computation 
time for a given step-size than the other two formulas. 
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As previously stated, the input data chosen for this study resulted in a 
trajectory that has both nominal and highly abnormal vehicle dynamics. The 
first 70 seconds of the trajectory represents normal dynamics since the 
vehicle is disturbed only by a wind. An engine failure in combination with the 
wind resulted in dynamics which very closely approached control loss following 
70 seconds. Thus, the formula 4-3 has proven to be the most practical formula 
for integrating both nominal and highly off-nominal vehicle dynamics. 

With the 4-3 formula established as a practical one to use in the liftoff 
program, an investigation of the value of the formula 4-3 regulator for control- 
ling the step-size in the program becomes a logical extension of this analysis. 
As previously stated, twelve distinct values for the regulator can be obtained 
at each integration step in the liftoff program. Effective control of the step- 
size requires a regulator that is sensitive to the total accuracy of the integra- 
tion process. Since all the parameters obtained by integration are related 
through the differential equations, the selected regulator must be a function 
of the regulators for the parameters that contain the largest percentages of 
integration error. Here, the percentage of integration error for a parameter 
is defined to be: 100 x integration error r absolute magnitude of the param- 
eter. By assuming a regulator to be a good approximation of the actual integra- 
tion error, we can determine these parameters by analyzing the twelve regula- 
tor values for various integration step-sizes. Figures 19 through 30 present 
the time variations of the twelve difference regulators for integration step- 
sizes from . 03125 to . 50 . The Saturn V dynamic data used to generate 
Tables 1, 2, and 3 were also used to generate these figures. 
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Table 8 lists the twelve parameters that are integrated and gives the 


largest percentage of error associated with each one. The largest percentage 
of error for a parameter was computed by using the peak value of the param- 
eter’s regulator at a step-size of . 50 as the maximum integration error. 

Table 8 indicates that the regulators for co , co , and c o are the ones most 

& x y z 

sensitive to the total accuracy of the integration process. A comparison of the 

regulators for these three parameters (Figures 25, 26, and 27) with the other 

nine regulators verifies this. That is, the response of these three regulators 

to the engine failure and wind gust disturbances at 70 seconds was much 

greater than the response of the other nine regulators. 

With the regulators for a; , co , and co established as the ones most 
& x y z 

sensitive to the total integration accuracy, the question of how to combine 

them into one total regulator must now be answered. Figures 25, 26, and 27 

show that the regulator for co^ is approximately ten times larger than the 

regulators for co^ and co^ . This prevents the three regulators from being 

combined directly. Equation (3-33) shows that any disturbance that affects 

the roll attitude also affects the pitch and yaw attitude. (Notice that 

co & &<p /dt . co ~ dcp /dt , and co ~ da? /dt .) Based on these two facts, 
x Y r‘ y P z y 

the decision was made to use the root-sum -square of the regulator for co^ 
and co as the total regulator for controlling the step-size in the liftoff pro- 
gram. Figure 31 presents the time variation of the total regulator for step- 
sizes from . 03125 to .50. 
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Limits must be established for this regulator in order to use it to 


control the integration step-size. These limits are used in combination with 
maximum and minimum restrictions on the step-size to maintain an acceptable 
integration error. The step-sizes required to obtain an acceptable solution 
can be determined from Tables 1, 2, and 3. From liftoff to 70 seconds, a 
step-size of . 50 is acceptable; from 70 to 90 seconds a step-size between 
. 125 and . 25 is required. Thus, the maximum and minimum restrictions on 
the step-size would be . 50 and . 125, respectively. The upper limit, L 2 , 
for the total regulator can be established from Figure 31. Before 70 seconds, 
the regulator for a step-size of . 50 does not exceed . 00025. Following 70 
seconds, the regulator exceeds . 00025 once for a step-size of . 125 and twice 
for a step-size of . 25. Therefore, L 2 was chosen to be . 00025. An effective 
value for the lower limit, , was established from a tabulation of the time 
variation of the total regulator with step-size. A value of . 00005 was chosen 
for hi since this value would allow the step-size to increase only when the 
regulator began to reach a significant peak. These values for L-^ and L 2 
were verified as being reasonable by a number of computer runs using small 
variations to these numbers. The Saturn V dynamic data used to generate 
Tables 1, 2, and 3 were also used in these verification runs. Since the above 
regulator and step-size limits were established for a trajectory that contains 
both nominal and highly disturbed dynamics, they should be valid for any 
trajectory required by a vehicle design or design assurance study. 
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The formula 4-3 routine for the liftoff program was modified to use 
the total regulator, R , for controlling the step -size. If R < L* , the 
step-size is doubled; if R > L 2 , the step-size is halved; if Lj < R < L 2 , 
the step-size is not changed. However, before the step-size is changed, the 
value to which the step-size is to be changed is checked to make sure it is not 
outside the step -size restrictions. 

With Li = „ 00005 , L 2 = . 00025 and . 50 and . 125 as the maximum 
and minimum step-size, respectively, a trajectory was computed with the 
modified formula 4-3 routine. The solution of this trajectory at the 70, 80, 
and 90 second time points is presented by Table 9. This solution is not only 
more accurate than the formula 4-3 solution with a step-size of .25 (Tables 
1, 2, and 3) , it also required 25 percent less computer time. Thus, the 
chosen regulator has been shown to be an effective indicator for controlling 
the integration error and reducing the computer time. 

This chapter has established that the fourth-order formula 4-3 is a 
practical integration formula to use in the liftoff program. Since this program 
contains the basic differential equations of motion of which most rigid-body 
trajectory programs are composed, the formula 4-3 numerical integration 
routine has been shown to be more practical to use in these programs than 
either the fourth-order Runge-Kutta formula or the fifth-order formula 5-5. 

In addition, the technique for controlling the integration step-size associated 
with the formula 4-3 has been shown to effectively increase the efficiency of 
the formula by as much as 25 percent when the trajectory being computed 


contains large dynamic disturbances. 


81 



00 

N> 


Table 1. Solution of Trajectory at 70 Seconds 


Step -Size 

<p 

p 

<P 

y 

cp 

r r 

xs 

p 

yS p 

zs 

P 

Time 51 ' (min) 

Formula 4-3 

. 03125 

-28.7873 

-.4728 

. 0027 

8947. 70 

11009. 18 

36465. 81 

23.646 

. 0625 

-28. 7873 

-.4728 

.0027 

8947.69 

11009. 17 

36465. 80 

12.066 

. 125 

-28. 7873 

-.4725 

. 0027 

8947. 18 

11009.17 

36465. 73 

6.306 

.25 

-28. 7873 

-.4729 

. 0027 

8947. 88 

11009.14 

36465. 72 

3.321 

. 50 

i 1 

-28.7872 

-.4717 

.0027 

8946.41 

11009. 13 

36465. 54 

2. 005 

Runge-Kutta Formula 

.03125 

-28.7873 

-.4726 

.0027 

8947.34 

11009.19 

36465. 79 

30.814 

.0625 

-28. 7873 

-.4726 

.0027 

8947.38 

11009. 18 

36465. 78 

15.650 

. 125 

-28. 7873 

-.4724 

.0027 

8947. 15 

11009. 18 

36465. 72 

8.075 

.25 

-28. 7874 

-.4727 

.0027 

8948. 00 

11009. 14 

36465.69 

4.258 

.50 

-28.7873 

-. 4709 

.0027 

8946.41 

11009. 14 

36465. 43 

2.432 


Formula 5-5 


. 03125 

-28. 7876 

-.4743 

.0027 

8950. 43 

11009. 00 

36465.65 

42.697 

.0625 

-28.7876 

-.4741 

.0027 

8950. 10 

11009.02 

36465.75 

21.431 

. 125 

-28. 7876 

-.4740 

. 0027 

8950. 11 

11009.06 

36465.75 

10. 820 

.25 

-28. 7875 

-.4732 

.0027 

8949. 16 

11009.05 

36465.65 

5.663 

.50 

-28. 7873 

-.4719 

. 0027 

8947. 76 

11008. 96 

36465. 74 

3. 185 


* Actual computer time used for computing the first 90 seconds of the trajectory. This number does not 
include time used in loading deck or reading data. 



Table 2. Solution of Trajectory at 80 Seconds 


Step -Size 

p 

cp 

y 

T r 

xs 

p 

y s 

P 

zs 

P 

Time* (min) 

Formula 4-3 

.03125 

-38.7233 

-16.5128 

-12.4525 

12180.78 

12036.13 

42546.21 

23.646 

.0625 

-38.7237 

-16.5187 

-12.4552 

12180.74 

12036. 11 

42546.19 

12.066 

.125 

-38.7233 

-16.5169 

-12.4545 

12180. 19 

12036. 18 

42546.11 

6.306 

.25 

-38.7158 

-16.4344 

-12.4212 

12181.15 

12036.22 

42546. 14 

3.321 

.50 

-39.0284 

-18.5710 

-12.8615 

12177.08 

12031.53 

42545.32 

2.005 

Runge-Kutta Formula 

.03125 

-38.7211 

-16.5001 

-12.4479 

12180.36 

12036. 18 

42546. 18 

30.814 

.0625 

-38.7266 

-16.5366 

-12.4628 

12180. 38 

12036.11 

42546. 15 

15.650 

.125 

-38.7215 

-16.4964 

-*12 . 4460 

12180. 17 

12036.24 

42546. 10 

8.075 

.25 

-38.7129 

-16.3931 

-12.4018 

12181. 32 

12036.36 

42546.13 

4.258 

.50 

-38. 9201 

-17.8007 

-12.7405 

12177.60 

12032. 99 

42545.28 

2.432 

Formula 5-5 

.03125 

-38.6932 

-16.2973 

-12.3574 

12184.00 

12036.10 

42546.09 

42.697 

.0625 

-38.6962 

-16.3158 

-12.3661 

12183.63 

12036.12 

42546.23 

21.431 

.125 

-38.6922 

-16.2869 

-12.3540 

12183.64 

12036.23 

42546.25 

10.820 

.25 

-38.6869 

-16.2431 

-12.3358 

12182.59 

12036.48 

42546. 14 

5.663 

.50 

-38.6550 

-15. 9663 

-12. 1963 

12181.28 

12037.31 

42546.33 

3. 185 


* Actual computer time used for computing the first 90 seconds of the trajectory. This number does not 
include the time used in loading the deck or reading data. 
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Table 3. Solution of Trajectory at 90 Seconds 


Step -Size 

<p 

p 


(p 

r 

xs 

p 

yS p 

zs 

p 

Time* (min) 

Formula 4-3 

.03125 

-44. 9268 

15. 9817 

7.6516 

15575.30 

12678.83 

49598.65 

23.646 

.0625 

-44. 9276 

16. 0061 

7.6691 

15575.20 

12678. 51 

49598.58 

12. 066 

. 125 

-44. 9273 

15. 9925 

7.6596 

15574.61 

12678.68 

49598.49 

6.306 

.25 

-44. 9169 

15.6544 

7.4057 

15576. 48 

12682. 91 

49599.23 

3.321 

. 50 

-44.3560 

19.4485 

9. 7651 

15545. 42 

12518.26 

49581. 74 

2.005 


Runge-Kutta Formula 


.03125 

-44.9252 

15. 9302 

7.6136 

15574.95 

12679. 55 

49598.69' 

30. 814 

.0625 

-44. 9300 

16.0834 

7. 7227 

15574.63 

12677.58 

49598.39 

15.650 

. 125 

-44. 9254 

15.9176 

7.6010 

15574. 76 

12679. 82 

49598. 62 

8.075 

.25 

-44.9133 

15. 5145 

7.2618 

15577.00 

12685. 16 

49599. 52 

4.258 

.50 

-44.8984 

20.0205 

9. 9560 

15557. 96 

12592.42 

49587. 85 

2. 432 


Formula 5-5 


.03125 

-44. 8996 

15.0857 

6. 8810 

15580.62 

12689.61 

49600. 12 

42.697 

.0625 

-44. 9021 

15. 1658 

6. 9562 

15580.09 

12688.74 

49600. 16 

21.431 

. 125 

-44. 8990 

15.0514 

6. 8434 

15580.31 

12690.27 

49600. 40 

10.820 

.25 

-44. 8951 

14. 8850 

6.6719 

15579.45 

12692.71 

49600. 52 

5.663 

. 50 

-44. 8646 

13. 8237 

5. 8255 

15580. 10 

12706.29 

49602.60 

3.185 


* Actual computer time used for computing the first 90 seconds of the trajectory. This number does not 
include the time used in loading the deck or reading data. 



Table 4. One -Sigma Data Tolerances 


Tolerance Number Tolerance 


1 . 1 degree thrust vector misalignment of each 
S-IC stage engine in the same direction 

2 . 3 percent variation of the total thrust magnitude 
of the S-IC stage engines 

3 2/3 inch lateral deviation of the predicted 
location of the missile center-of-gravity 

4 2/3 meter variation of the predicted aerodynamic 
center -of-pressure 


5 

6 
7 


2 percent variation of the predicted aerodynamic 

normal force coefficient C _ 

N 

3. 33 percent variation of the predicted 
aerodynamic drag coefficient 

2722 kilogram variation in the predicted initial 
mass of the missile 
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Table 5. Effects of One -Sigma Data Tolerances on Solution of Trajectory at 70 Seconds 


Tolerance* 

Number 

cp 

p 

y 

<p 
^ r 

xs 

p 

YS P 

zs 

p 

+ 1 

-28.7017 

-.3893 

.0031 

8949. 35 

11054.09 

36401.81 

-1 

-28. 8730 

-.5559 

.0022 

8944. 70 

10964. 07 

36529.69 

+2 

-28. 8044 

-. 5687 

.0029 

9062.40 

11007.56 

36479. 13 

-2 

-28.7734 

-.3990 

.0033 

8831.88 

11010.71 

36452.34 

+3 

-28. 8178 

-.4992 

-.0347 

8945.35 

10974. 86 

36514.07 

-3 

-28.7569 

-.4461 

.0400 

8948. 84 

11043.34 

36417.43 

+4 

-28.7886 

-.4015 

.0032 

8947. 15 

11013.52 

36465.78 

-4 

-28. 7860 

-. 5442 

. 0022 

8947.22 

11004. 87 

36465.68 

+5 

-28. 7874 

-.4647 

.0027 

8947. 17 

11010. 91 

36465.73 

-5 

-28. 7871 

-.4803 

.0028 

8947.20 

11007.46 

36465.73 

+6 

-28.7895 

-. 4849 

. 0027 

8962.74 

11008. 95 

36467. 89 

-6 

-28.7852 

-. 4613 

. 0028 

8933. 03 

11009.34 

36463.56 

+7 

-28. 7925 

-.4997 

.0026 

8983. 45 

11008.58 

36470.28 

-7 

-28.7824 

-.4473 

. 0029 

8910. 99 

11009. 76 

36461.20 


* The tolerances are defined by Table 4. A positive tolerance number indicates that the tolerance w.as applied 
in the most helpful way; a negative tolerance number indicates that the tolerance was applied in the most 
harmful way. 



Table 6. Effects of One-Sigma Data Tolerances on Solution of Trajectory at 80 Seconds 


Tolerance* 

Number 

cp 

P 


(f 

r 

xs 

p 

yS p 

zs 

p 

+i 

-38.3341 

-14. 9379 

-11.5715 

12184. 76 

12094.26 

42468.09 

-1 

-39.1924 

-18.3642 

-12.5339 

12175.13 

11977.75 

42624.17 

+2 

-38.4141 

- 9.7470 

- 6.2857 

12334.38 

12040.72 

42573.47 

-2 

-40.1522 

-26.3994 

-13.3281 

12025.63 

12032.69 

42518.46 

+3 

-38.6592 

-16. 1072 

-12.4730 

12177.77 

11994. 18 

42604.66 

-3 

-38.8017 

-16.9819 

-12.3504 

12182.35 

12077.95 

42487.61 

+4 

-38.1799 

- 9.4252 

- 5. 8506 

12185.03 

12057. 55 

42546. 86 

-4 

-41 . 5404 

-32.0683 

-13.3339 

12172.33 

12010.78 

42544.29 

+5 

-38.3416 

-13. 8243 

-10.6210 

12181. 95 

12043.89 

42546. 19 

-5 

-39.3255 

-20.0300 

-12.8337 

12178.20 

12028. 17 

42545. 95 

+6 

-38. 5942 

-15.3986 

-11. 9084 

12203. 54 

12036.76 

42551.62 

-6 

-38.7556 

-16. 9253 

-12.6127 

12159. 15 

12037. 06 

42540. 83 

+7 

-38.4446 

-14.0336 

-10.8898 

12229.68 

12037.04 

42555.64 

-7 

-39.0920 

-19.2638 

-12.8763 

12130. 81 

12035.54 

42536.59 


* The tolerances are defined by Table 4. A positive tolerance number indicates that the tolerance was applied 
in the most helpful way; a negative tolerance number indicates that the tolerance was applied in the most 
harmful way. 
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Table 7. Effects of One-Sigma Data Tolerances on Solution of Trajectory at 90 Seconds 


Tolerance v 
Number 

% 

<p 

y 

(P 

r 

xs 

P 

yS p 

zs 

P 

+ 1 

-44. 7510 

10.3892 

4.7007 

15592. 16 

12816.49 

49514.55 

: -l 

-44.2339 

18.2383 

9. 0583 

15543.26 

12455.62' 

49678.34 

: +2 

-44. 1665 

2.7922 

. 9414 

15798.29 

12862. 91 

49682.44 

-2 

-42. 7376 

-37.5092 

-407. 550 

14636. 99 

12306.83 

48931.55 

: +3 

-44. 9339 

14. 5655 

6. 7857 

15574. 91 

12647.18 

49699. 56 

-3 

-44. 9066 

17.6410 

8. 1382 

15573. 06 

12704.55 

49526. 97 ' 

+4 

-44. 1551 

1. 9187 

.5546 

15614.33 

12890. 09 

49635. 84 

-4 

-50.6667 

-46.7982 

-466. Oil 

14696. 97 

12333. 00 

48871.73 

+5 

-44.6586 

8. 0330 

3.5145 

15593.61 

12787.78 

49614. 93 

-5 

-209. 788 

-17. 8573 

-289. 193 

15343.07 

12191. 16 

49450. 16 

+6 

-44. 8128 

11. 8636 

5.6754 

15614. 90 

12727.34 

49616.93 

-6 

-44. 9548 

17.4259 

8.4833 

15541.61 

12659.01 

49584. 86 

+7 

-44. 6642 

8. 8958 

3. 9456 

15652.43 

12771. 53 

49630.58 

-7 

-42.6074 

10.6717 

5.3222 

15463. 93 

12418.64 

49556. 72 


* The tolerances are defined by Table 4. A positive tolerance number indicates that the tolerance was applied 
in the most helpful way; a negative tolerance number indicates that the tolerance was applied in the most 
harmful way. 



Table 8, Integration Error Percentages 


Maximum Percentage 
of 

Parameter Integration Error 


XSp 

.38 x 

10"' 

dX /dt 
s 

. 12 x 

10° 

YSp 

.29 x 

io -: 

dY /dt 
s 

. 22 x 

10° 

ZSp 

. 12 x 

10 -: 

dZ /dt 
s 

X 

00 

(N 

10° 

CO 

X 

CO 

X 

10 1 

CO 

y 

. 19 x 

10 1 

CO 

z 

. 20 x 

10 1 

cp 

p 

. 93 x 

i<r : 

<p 

y 

. 77 x 

10° 


. 82 x 

10° 
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Table 9. Formula 4-3 Solution of Trajectory with Variable Step-Size 


Trajectory 
Time (sec) 



<P 

r 

XS 

P 

YS 

P 

ZS 

P 

Computer 
Time (min) 

70 

-28.7874 

- .4726 

.0027 

8947. 90 

11009. 12 

36465.63 


80 

-38. 7193 

-16. 4557 

-12.4490 

12181.17 

12036. 17 

42546.04 


90 

-44.9181 

15. 7993 

7.5117 

15576.26 

12681.74 

49598. 98 

2.466* 


* Actual computer time used for computing the first 90 seconds of the trajectory. This number does not 
include time used in loading deck or reading data. 
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